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Abstract 

A wind tunnel experiment for characterizing the 
aerodynamic and propulsion forces and moments 
acting on a research model airplane is described. The 
model airplane, called the Free-flying Airplane for 
Sub-scale Experimental Research (FASER), is a 
modified off-the-shelf radio-controlled model 
airplane, with 7 ft wingspan, a tractor propeller 
driven by an electric motor, and aerobatic capability. 
FASER was tested in the NASA Langley 12- foot 
Low- Speed Wind Tunnel, using a combination of 
traditional sweeps and modem experiment design. 
Power level was included as an independent variable 
in the wind tunnel test, to allow' characterization of 
power effects on aerodynamic forces and moments. 
A modeling technique that employs multivariate 
orthogonal functions was used to develop accurate 
analytic models for the aerodynamic and propulsion 
force and moment coefficient dependencies from the 
wind tunnel data. Efficient methods for generating 
orthogonal modeling functions, expanding the 
orthogonal modeling functions in terms of ordinary 
polynomial functions, and analytical orthogonal 
blocking were developed and discussed. The 
resulting models comprise a set of smooth, 
differentiable functions for the non-dimensional 
aerodynamic force and moment coefficients in terms 
of ordinary' polynomials in the independent variables, 
suitable for nonlinear aircraft simulation. 
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Nomenclature 

a parameter vector 

C L ,Cd’ Cy lift, drag, and side force coefficients 

Q rolling moment coefficient 

C m pitching moment coefficient 

C„ yawing moment coefficient 

Cov covariance matrix 

J cost function 

MDOE Modem Design Of Experiments 

n number of model terms 

N total number of data points 

OF AT One Factor At a Time 

PSE predicted squared error 

pwr power level, percent 

p modeling function vector 

T thrust force, lbf 

x independent variable vector 

y measured output vector 

a angle of attack, deg 

/3 sideslip angle, deg 

S Q aileron deflection, deg 

S e elevator deflection, deg 

Sf flap deflection, deg 

S r rudder deflection, deg 

o 2 variance 

ordinary' polynomial function vector 

superscripts 

T transpose 

estimate 

-1 matrix inverse 

normalized 
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subscripts 


max 

maximum 

min 

minimum 

o 

nominal 


Introduction 

Modem aeronautical research involves expanded 
flight envelopes, as a result of the desire for improved 
fighter maneuverability for tactical advantages, and the 
desire to improve flight safety. The expanded flight 
envelopes involve nonlinear aerodynamics which must 
be modeled accurately. 

Since nonlinear aerodynamics are much more 
complex than linear aerodynamics, more sophisticated 
experimentation is required to accurately characterize 
the functional dependencies. Nonlinear aerodynamics 
violate linear modeling assumptions such as 
superposition, quasi-steady flow, and no 
interdependence of the effects of states and controls. In 
addition, aircraft designs have evolved with increasing 
numbers of control effectors. Traditional wind tunnel 
testing methods would set each control effector to 
different fixed levels while sweeping through angle of 
attack and sideslip angle, for example. With such an 
approach, the number of data points required increases 
exponentially with the number of control effectors, if 
information on control surface interaction effects is 
desired. These considerations highlight the need to 
develop more efficient wind tunnel testing and 
modeling techniques to accurately characterize 
nonlinear aerodynamics, with possible interaction 
effects among a large number of independent variables. 

At NASA Langley, the Free-flying Airplane for 
Sub-scale Experimental Research (FASER) is being 
developed to study problems such as that described 
above. FASER is a modified off-the-shelf radio control 
model called the Ultra-Stick, manufactured by Hangar 
9, Ltd., see Figure 1. FASER has a conventional high- 
wing and tail configuration with 7 ft wingspan, a 
foldable tractor propeller driven by an electric motor, 
and aerobatic capability. 

The purpose of FASER is to provide an 
inexpensive aircraft for developing and demonstrating 
advanced experiment design, data analysis and 
modeling techniques, and control law' design methods. 
As long as the goal is technology demonstration or 
basic research, a sub-scale model that is not 
dynamically scaled for a specific full-scale aircraft is a 
completely acceptable test vehicle for these purposes. 


FASER was designed so that the flight vehicle 
could be installed in the wind tunnel, see Figure 1. 
This avoids any Reynolds number or scaling effects, 
and ensures that there are no physical differences 
between the wind tunnel model and the flight vehicle. 

In contrast, full scale flight tests and drop model 
tests are expensive and are sometimes separated by 
months (and even years) for a particular research 
activity. The number of these tests is always tightly 
constrained by budget. There is also a substantial 
difference in the cost of overhead if the aircraft is to be 
kept in flyable condition. Since FASER is inexpensive 
and unmanned, risks can be taken in research and 
development that could never be tolerated in a piloted 
flight test or even in a drop model test. Advances in 
instrumentation now make it possible to instrument a 
sub-scale model aircraft with research-quality, 
miniaturized flight test instrumentation at a reasonable 
cost. 

This paper describes the experiment design, data 
analysis, and mathematical modeling involved in 
developing a wind tunnel database for the FASER 
aircraft. Accuracy of this database is critical for the 
development of a high-fidelity nonlinear simulation to 
be used for control law design, flight envelope 
expansion, flight experiment design, and pilot training. 
A preliminary version of the nonlinear simulation for 
FASER has already been developed, using U.S. Air 
Force DATCOM to generate an aerodynamic model, 
with experimentally-determined values for the mass 
and inertia characteristics of FASER 1 . The work 
described in this paper will upgrade the aerodynamic 
model with analytic models derived from wind tunnel 
data, add an engine thrust model, and include 
propulsion effects on the aerodynamics. Since FASER 
was intended to be a research vehicle from the outset, 
the approach used for the experiment design and data 
analysis for the wind tunnel testing was not traditional. 
This paper explains how the wind tunnel testing was 
done, and examines the results. The paper also 
describes a method for generating orthogonal modeling 
functions based on the independent variable data, along 
subsequent expansion of the orthogonal modeling 
functions in terms of ordinary' multivariate 

polynomials. This method is slightly different from 
that described in Refs. [2] and [3], and represents an 
evolutionary improvement of the technique. In Ref. 
[2], the concept of response surface modeling using 
multivariate orthogonal functions was successfully 
applied to inference subspaces for limited ranges of 
angle of attack and Mach number with fixed sideslip 
angle and control surface deflections. This paper 
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extends the multivariate orthogonal function modeling 
concept to identify aerodynamic models for a large 
flight envelope, with more independent variables. 

Experiment Design 

For this wind tunnel test, the fundamental 

objective was to find a mathematical description for the 

dependence of non-dimensional aerodynamic force and 
moment coefficients on independent variables that are 
varied during the experiment. Each mathematical 
description or model can be thought of geometrically as 
a hyper-surface, also called a response surface. Critical 
issues for successfully identifying an adequate response 
surface model from experimental data include the 

experiment design (or, how the independent variable 
values are set when measuring the output responses), 
noise level on the measured outputs, identification of a 
mathematical model structure that can capture the 

functional dependence of the output variables on the 
independent variables, accurate estimates of unknown 
parameters in the identified model structure, and the 
ability of the identified model to predict outputs for 
data that was not used to identify the response surface 
model. 

The experiment design used for FASER wind 
tunnel testing was a hybrid design broken down into a 
series of procedures. The procedures are listed in 
Table 1. Randomization 3 4 ' 6 was used throughout the 
testing, to separate independent variable effects from 
time-dependent systematic errors. 

The first procedure consisted of randomized 
engine power sweeps with the wind tunnel air off, to 
determine the static thrust from the electric motor and 
the propeller. All of these runs were made with the 
model at zero angle of attack and zero sideslip angle, so 
that the thrust measurement was obtained from the 
longitudinal force measured by the balance mounted in 
the model. Figure 2 shows the static thrust plotted as a 
function of pulses per second from a Hall effect 
transducer on the electric motor, which is proportional 
to the propeller RPM. The model shown in Figure 2 is 
the result of a simple least squares fit of thrust to pulse 
count, using a quadratic model structure. This model 
structure was identified automatically from the data, 
using the orthogonal function modeling technique 
described later. 

In the second and third procedures, the approach 
was to use One Factor At a Time (OF AT) sweeps, 
wherein one independent variable is changed with all 
others held constant, to characterize the general 
topology of the response surfaces for the 


non-dimensional coefficients, get an idea of the 
response levels, collect basic static stability and trim 
information, and define the boundaries of the 
independent variable subspaces to be used for further 
experimentation. OF AT sweeps were used because the 
data acquisition system in the NASA Langley 12-foot 
Low-Speed Wind Tunnel is set up to collect OFAT data 
in an automated fashion, making the sweeps very 
efficient in terms of collecting data points in minimum 
time. However, the initial OFAT sweeps are really 
only intended for qualitative use, namely to define the 
boundaries for subspaces that will be the focus of 
detailed experimentation and modeling in procedure 
four. One advantage of operating in this manner is that 
any issues related to instrumentation, data collection, or 
experimental procedures can be worked out during the 
OFAT sweeps without impact on the experiment, 
because the data from the OFAT sweeps is being used 
for qualitative purposes only. This approach also 
provides a good rough overview^ of the landscape that 
characterizes the dependence of force and moment 
coefficients on the independent variables. 

The independent variables for the FASER wind 
tunnel test were angle of attack a, sideslip angle 
power level /ntr, elevator deflection S e , aileron 

deflection 5 a , rudder deflection S r , and flap 
deflection 8 f . The response variables were 

non-dimensional aerodynamic coefficients for lift, drag, 
and side forces (Q ,Q>, and C Y ) . and rolling, 

pdching, and yawing moments (Cj , C m , and C„ ) . 

Each data point produced measured values for all 
independent and response variables. 

The experiment was designed assuming any of the 
independent variables could influence any of the 
response variables. It was assumed (as an initial guess 
only) that the dependencies could be modeled 
accurately with polynomial terms in the independent 
variables of order 3 or less within each independent 
variable subspace. In addition, it was assumed that 
longitudinal controls ( S e ,Sf , and pwr) do not interact 

w ith lateral/directional controls ( S a and S r ). In 
practical terms, this meant that longitudinal and 
lateral/directional controls were not varied 
simultaneously to allow their mutual interaction effects 
to be quantified. As a result, the subspaces were called 
’'longitudinal” if the longitudinal controls were moved, 
and "lateral/directional” if the lateral/directional 
controls were moved. 

Based on experience with similar airplanes, it was 
known that the dependence of non-dimensional 
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aerodynamic force and moment coefficients on control 
surface deflections could be modeled with low order 
polynomials for the entire range of control surface 
deflections. With that assumption, it was not necessary 
to vary' the control surface deflections to search for 
inference subspace boundaries along the dimensions of 
the independent variable space corresponding to control 
surface deflections. Inference subspace boundaries 
were therefore sought only for a, /?, and pwr. These 
independent variables were varied using OF AT sweeps 
to identify the inference subspace boundaries. 

Figure 3 shows an OFAT sweep on angle of 
attack. The vertical lines mark the selected subspace 
boundaries in angle of attack, which are intended to 
mark the boundaries of regions where the character of 
the response surfaces change. There is a trade-off in 
selecting the subspace boundaries, in that more 
subspaces mean more individual experimentation 
regions in procedure four, while fewer subspaces 
generally require more resources in the data collection 
and modeling for each subspace. Figures 4 and 5 show 
OFAT sweeps on sideslip angle and power level, with 
the selected subspace boundaries marked by vertical 
lines. All control surface deflections were zero for the 
sweeps shown in Figures 3-5. 

The full independent variable ranges tested are 
listed in Table 2. Tables 3 and 4, and corresponding 
Figures 6, 7, and 8, show the definitions of the 
inference subspaces in terms of boundary values of 
angle of attack, sideslip angle, and power level. The 
symmetry of the vehicle was used to justify omitting 
testing in most of the subspaces with high negative 
sideslip angles, see Figures 6, 7, and 8. Ail the 
subspaces together comprised the full inference space, 
defined by the full range of the independent variables in 
Table 2. All control surfaces were tested over their full 
physical deflection ranges for each subspace. 

In the fourth and final procedure. Modem Design 
Of Experiment (MDOE) techniques 4 * ' 6 were applied to 
each defined subspace in order to obtain the most 
accurate and complete characterization of the functional 
dependencies, and also to make sure all interaction 
effects were adequately modeled. Refs. [2], [4]-[6] 
describe and demonstrate in detail the advantages of the 
MDOE approach compared to traditional OFAT for 
detailed modeling of the functional dependencies, both 
in terms of the modeling accuracy and in the economy 
of experimentation resources required to arrive at an 
acceptable result. 

Within each subspace, independent variables were 
set according to normalized values. Normalized 


independent variable values are found by mapping the 
independent variable values in engineering units for 
each subspace onto the interval [-1,1]. The 
normalization of each independent variable was 
implemented by 

x = -1 + 2 -S —— - ( 1 ) 

l X max ~ x min ) 

where x was the normalized value of the independent 
variable, and the independent variable range in 
engineering units was [ x min , x max ] . The inverse 
transformation was 

(x + l) 

x ~~ x min 4 2 \ X max ~~ X min ) ( 2 ) 

Ail modeling for the inference subspaces was 
done using normalized values of the independent 
variables. The final models used for prediction were 
written in terms of engineering units. 

Second-order central composite design 5 6 in four 
independent variables (for the power-off subspaces) or 
five independent variables (for the power-on 
subspaces) was used in each subspace, augmented with 
a 3 rd order D-optimal design 5 6 in the appropriate 
number of independent variables. A two-dimensional 
projection of this constellation of data points in 
normalized independent variable space is shown in 
Figure 9. The central composite design points occupy 
the comers, the centers of each face, and the center of 
the normalized subspace, while the D-optimal points 
generally fill in between. Although some of the 
D-optimal points are coincident with the central 
composite design points, the number of times that this 
happens is not represented accurately in Figure 9, 
because of the projection onto two dimensions. This 
experiment design enabled identification of up to 3 rd 
order functional dependencies and interaction effects. 
Provisions were made to augment the designed 
experiment if the data indicated a lack of fit that 
required modeling functions with higher than 3 rd order. 

Instrumentation and Data Collection 

FASER was used as the wind tunnel model and 
tested at a nominal flight speed, thus avoiding scaling, 
Reynolds number, or geometric dissimilarity issues for 
comparisons with future flight test data. All control 
surfaces were instrumented with potentiometers. Air 
data vanes and airspeed pinwheels were installed on 
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booms attached at each wing tip and extending 1 chord 
length in front of the wing. The air data sensors, which 
will be used for flight testing, were calibrated as part of 
the wind tunnel experiment, since aerodynamic 
incidence angles and airspeed were carefully controlled 
and measured in the wind tunnel. The wind tunnel 
balance was installed near the c.g. of the airplane in the 
space normally occupied by the accelerometer and rate 
gyro package during flight test operations. 

Control surface deflections and power level were 
automatically set to the values required by the 
experiment design via a serial port interface from a 
laptop computer in the control room to a 
servomechanism controller onboard on the airplane. 
The same onboard equipment will be used to command 
the control surfaces and power level during flight 
operations. Angle of attack and sideslip angle were set 
from the control room using servomechanisms driving 
the movable C-strut in the test section. The angle of 
attack and sideslip angles were set automatically during 
OFAT sweeps, but had to be set manually for the 
MDOE data points, using joystick controllers and a 
measurement feedback to the control room. Dynamic 
pressure in the tunnel was regulated to 2 psf by an 
automatic closed loop control on the wind tunnel fan 
motor speed. 

The experimental set-up was designed to 
accommodate the MDOE approach, which typically 
requires changes in more than one independent variable 
for successive data points. Since the control surface 
and power level settings were automated using the 
laptop computer, each data point required only manual 
setting of angle of attack and sideslip angle using the 
joysticks in the control room. Each data point was 
taken as the average of a ten-second dwell using a 
sampling rate of 100 Hz. 

The data for power-on subspaces was collected by 
interleaving power-on points with power-off points 
from other subspaces, in order to keep the engine 
temperature at low levels for extended testing periods. 
This was necessary to avoid damage to the electric 
motor. A regulated DC power supply was used to 
power the electric motor, so that batteries would not 
have to be swapped in and out. The power-on 
subspaces were limited to relatively low angle of attack 
and sideslip angle (see Figure 7), because of excessive 
vibration of the wind tunnel rig and model for powered 
runs at high angles of attack and/or high sideslip 
angles. 

The wind tunnel experiment described above was 
conducted over 4 weeks in May 2002, in the NASA 
Langley 12-foot Low-Speed Wind Tunnel. 


Modeling 

Typically, once the experimental data are 
collected, polynomials in the independent variables are 
used to model the functional dependence of the output 
variables on the independent variables, and the model 
parameters are estimated from the measured data using 
least squares linear regression 5 6 . The question of 
which polynomial terms should be included in the 
model for a given set of data, called model structure 
determination, gets more difficult with increases in the 
number of independent variables, increased ranges for 
the independent variables, or increased complexity of 
the underlying functional dependency . 

Various hypothesis testing techniques 6 7 can be 
used to identify an adequate model structure, but these 
methods are iterative and require the involvement of an 
experienced analyst. Neural networks using radial 
basis functions with subspace partitioning, or back 
propagation with layered and interconnected nonlinear 
activation functions, have also been applied to the 
response surface modeling problem 8 . For this type of 
approach, there is a loss of physical insight and a 
danger of over-fitting the data, because the model 
structures contain many parameters, typically with no 
mechanism for limiting the size of the model other than 
the judgment of the analyst. 

In this work, a nonlinear multivariate orthogonal 
modeling technique 2 ’ 3 was used to model response 
surfaces from wind tunnel data. The technique 
generates nonlinear orthogonal modeling functions 
from the independent variable data, and uses those 
modeling functions with a predicted squared error 
metric to determine appropriate model structure. The 
orthogonal functions are generated in a manner that 
allows them to be decomposed without ambiguity into 
an expansion of ordinary multivariate polynomials. 
This allows the identified orthogonal function model to 
be converted to a multivariate ordinary polynomial 
expansion in the independent variables, which provides 
physical insight into the identified functional 
dependencies. 

The next section briefly describes the multivariate 
orthogonal function modeling approach. In the Results 
section, the multivariate orthogonal function modeling 
method is applied to identify response surface models 
for non-dimensional aerodynamic force and moment 
coefficients for inference subspaces, based on data 
from the FASER wind tunnel test. 
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Multivariate Orthogonal Functions 

Assume an TV-dimensional vector of response 
variable values, y = [y lt y 2 Y - m °deled in terms 

of a linear combination of n modeling functions 
Pj, j = 1,2 Each Pj is an TV-dimensional vector 

which in general depends on the independent variables. 
Then. 


Cov(a) = E^[a - a){a ~ a ) 1 j = a 2 j P J (8) 

where E is the expectation operator, and the error 
2 

variance a can be estimated from the residuals, 

v-y-Pd (9) 


y — tfj p\ + #2 Pi + ■ ■ ■ + a n Pn + £ (3) 

The Oj, y = l,2,...,w are constant model parameters to 

be determined, and e denotes the modeling error vector. 
Eq. (3) represents the usual mathematical model used to 
fit a response surface to measured data from an 
experiment. We put aside for the moment the 
important questions of determining how candidate 
modeling functions Pj should be computed from the 

independent variables, as well as which candidate 
modeling functions should be included in Eq. (3), 
which implicitly determines n. Now define an Nxn 
matrix P, 


P = [Pl,Pl, .Pn\ ( 4 ) 

and let a = \a x ,a 2 >.. -,a„f . Eq. (3) can be written as a 

standard least squares regression problem, 

y=Pa+£ (5) 

where y is a vector of measured dependent variable 
values, P is a matrix whose columns contain modeling 
functions of the measured independent variables, and a 
is a vector of unknown parameters. The variable £ 
represents a vector of errors that are to be minimized in 
a least squares sense. The goal is to determine a that 
minimizes the least squares cost function 

J= (y-Pa) T (y-Pa)=e 1 e (6) 

The parameter vector estimate that minimizes this cost 

function is computed from 3 ’ 5 " 6 7 

i= [ pTp ]~ 1 pT y 


The estimated parameter covariance matrix is 


cx 2 = 


(N-n) 


{y-Pa) T ( y-Pd ) 


v r v 

(N-n) 


( 10 ) 


and n is the number of elements in parameter vector a . 
Parameter standard errors are computed as the square 
root of the diagonal elements of the Cov[a) matrix 

from Eq. (8), using a from Eq. ( 1 0). 

Estimated model output is 

y = Pa (11) 


For response surface modeling, the modeling 
functions (columns of P) are often chosen as 
polynomials in the measured independent variables. 
This approach corresponds to using the terms of a 
Taylor series expansion to approximate the functional 
dependence of the output response variable on the 
independent variables. 

If the modeling functions are instead multivariate 
orthogonal functions generated from the measured 
independent variable data, advantages accrue in the 
model structure determination for response surface 
modeling. After the model structure is determined 
using the multivariate orthogonal modeling functions, 
each retained modeling function can be decomposed 
into an expansion of ordinary polynomials in the 
independent variables. Combining like terms from this 
final step puts the final model in the form of a Taylor 
series expansion. It is this latter form of the model that 
provides the physical insight, particularly in the case of 
modeling non-dimensional aerodynamic force and 
moment coefficients. This is the reason that aircraft 
dynamics and control analyses are nearly always 
conducted with the assumption of this form for the 
dependence of the non-dimensional aerodynamic force 
and moment coefficients on independent variables such 
as angle of attack and sideslip angle. 

Ref. [3] describes a procedure for using the 
independent variable data to generate multivariate 
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orthogonal modeling functions Pj , which have the 
following important property: 


p = rc _1 


(17) 


p]pj= 0 , i*j , i,j = l,2,...,n (12) 

It is also possible to generate multivariate 
orthogonal functions by first generating ordinary 
multivariate polynomials in the independent variables, 
then orthogonalizing these functions using 
Gram-Schmidt orthogonalization. The process begins 
by choosing one of the ordinary multivariate 
polynomial functions as the first orthogonal function: 

A=fi ( 13 ) 

where <f 3 is the ordinary multivariate polynomial 

function chosen to be the first orthogonal function. 
Then to make each subsequent ordinary multivariate 
polynomial function orthogonal to the preceding 
orthogonal function(s), define the 7 th orthogonal 
function p } as: 

H 

Pj^Sj-ZtYkjPk 7 = 2 « (14) 

k= 1 

where £ y is the 7 th ordinary multivariate polynomial 
vector, and the y k/ are scalars determined from 




pUj 

PkPk 


* = 1 , 2 7- 1 

7 = 2 « 


(15) 


where n is the total number of ordinary multivariate 
polynomials used as raw material for generating the 
multivariate orthogonal functions. Eq. (15) results 

T 

from multiplying both sides of Eq. (13) by p k , 
invoking the mutual orthogonality of p k , k=\ ,..., 7 , 
and solving for y kj . It can be seen from 

Eqs. ( 13)-( 1 5) that each orthogonal function can be 
expressed in terms of ordinary polynomial functions. If 
the pj vectors and the vectors are arranged as 

columns of matrices P and 5, respectively, and the 
y kj are elements in the k row r and 7 th column of an 

upper triangular matrix G with ones on the diagonal, 
then 


£ = PG 


(16) 


which leads to 


The columns of G ” 1 contain the coefficients for 
expansion of each column of P (i.e., each multivariate 
orthogonal function) in terms of the ordinary 
polynomial functions contained in the columns of £ . 
Eq. (17) can be used to express each multivariate 
orthogonal function in terms of ordinary multivariate 
polynomials. 

The orthogonal functions are generated in a 
manner that allows them to be decomposed without 
ambiguity into an expansion of ordinary multivariate 
polynomials. The orthogonalization process can be 
repeated using arbitrary ordinary multivariate 
polynomials to generate orthogonal functions of 
arbitrary’ order in the independent variables, subject 
only to limitations related to the information contained 
in the data. For the FASER wind tunnel data response 
surface modeling, the orthogonal modeling functions 
were generated in the manner described above. 

If an additional independent variable is introduced 
to represent blocking in the experiment, the orthogonal 
function modeling can be used to orthogonalize the 
block effects with respect to the other independent 
variable effects. A blocking variable is typically used 
to indicate some change in the experimental conditions 
that cannot be controlled by the experimenter. The 
blocking variable to the first power can simply be made 
one of the ordinary polynomial vectors . The 

orthogonalization procedure in Eqs. (13)-(15) makes 
the blocking variable orthogonal to the other 
orthogonal modeling functions, which are generated 
from ordinary multivariate polynomial functions. This 
approach allows arbitrary blocking of the data points in 
the experiment using analytical means alone. 

Normally, experiment designers try to arrange the 
normalized settings of the independent variables so that 
the blocking variable is orthogonal to the other 
independent variables and their polynomial 
combinations. This separates the block effect from the 
other model terms and allows identification of a block 
effect independent of the other terms in the model. 
However, with the analytical orthogonal blocking 
described above, the blocking variable is made 
orthogonal to the other modeling functions analytically, 
so that arranging the experiment so that independent 
variable settings and their polynomial combinations are 
orthogonal to the blocking variable is not necessary. 
All the models identified in this work included an 
analytic blocking variable that accounted for drift in 
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experimental conditions between the time when the 
central composite design points were collected and the 
time when the D-optimal design points were collected. 

With the modeling functions orthogonal, using 
Eqs. (4) and (12) in Eq. (7), the 7 th element of the 
estimated parameter vector a is given by 

°j =(pU)/(pj Pj) < 18 > 

Combining Eqs. (4), ( 6 ), and ( 1 1 )-( 1 2), 

j-r T y-i$(ry») < k » 

7=1 

or, using Eq. (18), 

= (20) 

7=1 

Eq. (20) shows that when the modeling functions 
are orthogonal, the reduction in the estimated cost 
resulting from including the term Oj Pj in the model 

depends only the dependent variable data y and the 
added orthogonal modeling function Pj . This 

decouples the least squares modeling problem, and 
makes it possible to evaluate each orthogonal modeling 
function in terms of its ability to reduce the least 
squares model fit to the data, regardless of which other 
orthogonal modeling functions are present in the 
model. When the modeling functions p } are instead 

polynomials in the independent variables (or any other 
non-orthogonal function set), the least squares problem 
is not decoupled, and iterative analysis is required to 
find the subset of modeling functions for an adequate 
model structure. 


*■> 

where <7 max is the maximum variance of elements in 
the error vector e , assuming the correct model 
structure, and n is the number of model terms. The 
PSE in Eq. (22) depends on the mean squared fit error 
j/N , and a term proportional to the number of terms 
in the model, n . The latter term prevents over-fitting 
the data with too many model terms, which is 
detrimental to model prediction accuracy 9 . The factor 
of 2 in the model over-fit penalty term accounts for the 
fact that the PSE is being used when the model 
structure is not correct, i.e., during the model structure 
determination stage. Ref. [9] contains further justifying 
statistical arguments and analysis for the form of PSE 
in Eqs. (21)-(22). Note that while the mean squared fit 
error J/n must decrease with the addition of each 
orthogonal modeling function to the model (by Eq. (19) 
or ( 20 )), the over-fit penalty’ term <7^ n/N increases 

with each added model term ( n increases). Introducing 
the orthogonal modeling functions into the model in 
order of most effective to least effective in reducing the 

mean squared fit error (quantified by a~ ( p^ p } j for 

the 7 th orthogonal modeling function) means that the 
PSE metric will always have a single global minimum 
value. Figure 10 depicts this graphically, using actual 
modeling results from Ref. [2]. Ref. [9] contains 
details on the statistical properties of the PSE metric, 
including justification for its use in modeling problems. 

For wind tunnel testing, repeated runs at the same 
test conditions are often available. If ol is the output 
variance estimated from measurements of the output for 
repeated runs at the same test conditions, then a^ ax 
can be estimated as 

= (23) 


The orthogonal modeling functions to be included 
in the model are chosen to minimize predicted squared 
error PSE, defined by 8 9 


pt;E Jy- p t) T (y-rt) ,„ t 2 « rn 

PSE + 2 <J max — (21) 


or 


PSE = — + 2al ax — 
N N 


( 22 ) 


If the output errors were Gaussian, Eq. (23) would 
correspond to conservatively placing the maximum 
output variance at 25 times the estimated value 
(corresponding to a 5a 0 maximum deviation). 

However, the estimate of may not be very good, 

because of relatively few repeated runs available or 
errors in the independent variable settings or drift errors 
when duplicating test conditions for the repeat runs. 
The 5<j 0 value was found to give the most accurate 

models in model identification algorithm testing done 
using simulated data. In Ref. [2], the model structure 
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determined using PSE was found to be virtually the 
same for a^ ulx in the range: 

9 <7~ < a 2 max < 100 ol (24) 

This happens because the plot of mean square fit 
error versus added orthogonal modeling functions is 
typically very flat in the region of minimum PSE , see 
Figure 10. 

Using orthogonal functions to model the response 
variable made it possible to evaluate the merit of 
including each modeling function individually as part 
of the model, using the predicted squared error PSE. 
Since the goal is to select a model structure with 
minimum PSE , and the PSE always has a single global 
minimum for orthogonal modeling functions, the model 
structure determination was a well-defined and 
straightforward process that could be (and was) 
automated. 

After the orthogonal modeling functions that 
minimized PSE were selected, each retained orthogonal 
function was expanded into an ordinary polynomial 
expression, and common terms in the ordinary 
polynomials were combined using double precision 
arithmetic to arrive finally at a multivariate model using 
only ordinary polynomials in the independent variables. 
Ordinary polynomial terms that contributed less than 
0.1 percent of the final model root-mean- square 
magnitude were dropped. 

Orthogonal modeling functions are useful in 
determining the model structures for the response 
variables using the PSE metric, by virtue of the 
properties of orthogonal modeling functions and the 
resultant decoupling of the associated least squares 
problem. The subsequent decomposition of the 
retained orthogonal functions is done to express the 
results in physically meaningful terms and to allow 
analytic differentiation for partial derivatives of the 
response variable with respect to the independent 
variables. 


Results 

Figures 11, 12, and 13 show' results for the 
response surface model fits to measured lift, drag, and 
pitching moment coefficient data obtained by applying 
the multivariate orthogonal function nonlinear 
modeling technique to experimental data from 
longitudinal subspace 16 for the FASER wind tunnel 
test. The crosses shown in the top plot of the figures 
arc measured data, and the circles are values computed 
from the identified response surface models. Values 

for al were found from twelve repeated runs at the 

normalized center point of the independent variable 
ranges, using the method described above. Model 
structure determination and parameter estimation were 
done automatically using the orthogonal function 
modeling technique. The orthogonal function modeling 
software allows manual override by the analyst in the 
model structure determination stage, but this option was 
not used for the results presented here. All data 
analysis and modeling was done on a 1.2 GHz personal 
computer running MATLAB® 6.5. The orthogonal 
modeling technique described above was implemented 
as an m-file function. 

The model residuals shown in the middle plot of 
Figures 11, 12, and 13 show a random character with 
magnitudes on the order of the noise level estimated 
from repeated runs. The dashed lines in the middle and 
lower plots in the figures represent ±cx 0 , the square 

root of the estimated output variance, computed from 
repeated runs at the subspace center point. The lower 
plots in the figures show that the prediction residual 
magnitudes are also within the estimated noise levels 
for the output responses. Results shown for this 
subspace were typical. 

Table 5 contains the identified model structures 
for C D , C L , and C m for this inference subspace. 

Based on the data, the identification algorithm 
determined that no sideslip angle dependence was 
needed for the longitudinal coefficients models in 
subspace 16, so the final models did not include 
sideslip angle. 

The results shown in Figures 11, 12, and 13 and 
the associated analysis of residuals indicated that the 
functional dependencies were accurately captured by 
the response surface models identified using the 
experiment design and modeling techniques described. 
Similar results were generated for the other longitudinal 
and lateral/directional subspaces listed in Tables 3 and 
4 and depicted in Figures 6, 7, and 8. 
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Figure 14 shows a screen capture of a 3D plotting 
tool developed to inspect the response surface fits to 
measured data. The symbols represent the measured 
data points, and the smooth surface is the identified 
response surface. The viewpoint for the 3D plot can be 
rotated and zoomed in or out, under the control of the 
analyst, and the two independent variables for the 3D 
plotting can be selected from the complete set of 
independent variables. This tool gives a good overview 
of the response surface topology, using 3D slices 
through the inference space. 

Concluding Remarks 

This paper describes and demonstrates an efficient 
and effective approach to developing a wind tunnel 
database for a research model airplane. The use of 
OFAT sweeps for subspace definition proved useful in 
defining subspace boundaries and as trial runs to 
identify and resolve problems related to experimental 
procedure, data collection, and instrumentation. The 
use of multivariate orthogonal modeling functions 
simplified the model structure determination problem 
and allowed arbitrary orthogonal blocking. The quality 
of the modeling and prediction results from the 
experiment design and modeling approach described 
and used were excellent. However, because the wind 
tunnel was not set up for MDOE experimentation, the 
overall efficiency of the testing was not what it could 
have been. 

Significant efficiency improvement would accrue 
if the data could be transferred to a computer for 
analysis in real time, and if the wind tunnel rig could be 
set automatically to arbitrary angles of attack and 
sideslip angles. In this case, the modeling software 
could automatically determine the model structure and 
parameter values in real time as each data point is 
taken, then compute model fit metrics and make 3D 
plots. The augmentation of the MDOE experiment 
designs for higher order model fits, or the re-definition 
of subspace boundaries could be automated as well, 
resulting in a testing operation where the human analyst 
would be needed only to provide high-level oversight 
of the entire operation. 
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Table 1 FASER Wind Tunnel Experiment Procedures 


Procedure 

Description 

1 

Power sweep with wind off for 
static thrust modeling 

2 

Angle of attack sweeps at min, max, and 
zero elevator, for static stability and trim 

3 

Angle of attack, sideslip angle, and power 
sweeps for inference subspace 
identification 

4 

2 nd order central composite design and 
3 rd order D-optimal design for inference 
subspace modeling 


Table 2 Independent Variable Ranges 


Variable 

a 

P 

pwr 

Se 

S a 

S r 

< 5 / 

Min 

-7.5 

-30 

0 

-25 

-25 

-30 

0 

Max 

80 

30 

100 

25 

25 

30 

30 

Units 

deg 

deg 

% 

deg 

deg 

deg 

deg 


Table 3 Longitudinal Inference Subspaces 


Inference 

Subspace 

Angle of 
Attack 
(deg) 

Power 

Level 

(percent) 

Sideslip 
Angle (deg) 


min 

max 

min 

max 

min 

max 

13 

-7.5 

10 

0 

0 

-10 

10 

23 

-7.5 

10 

0 

100 

-10 

10 

6 

-7.5 

10 

0 

0 

10 

30 

15 

10 

20 

0 

0 

-10 

10 

17 

10 

20 

0 

0 

10 

30 

24 

10 

20 

0 

100 

-10 

10 

16 

20 

40 

0 

0 

-10 

10 

18 

20 

40 

0 

0 

10 

30 

5 

40 

80 

0 

0 

-10 

10 

4 

40 

80 

0 

0 

10 

30 


Table 4 Lateral/Directional Inference Subspaces 


Inference 

Subspace 

Angle of 
Attack 
(deg) 

Power 

Level 

(percent) 

Sideslip 
Angle (deg) 


min 

max 

min 

max 

min 

max 

7 

-7.5 

10 

0 

0 

-10 

10 

8 

-7.5 

10 

0 

0 

10 

30 

26 

-7.5 

10 

0 

0 

-30 

-10 

19 

10 

20 

0 

0 

-10 

10 

20 

20 

40 

0 

0 

-10 

10 

21 

10 

20 

0 

0 

10 

30 

22 

20 

40 

0 

0 

10 

30 

11 

40 

80 

0 

0 

-10 

10 

10 

40 

80 

0 

0 

10 

30 


Table 5 Longitudinal Subspace Identified Model 


Inference 

Subspace 

Model Structure 

16 

C D = a o + a 2 + a 2 aS e + a 3 8f 

Ci = ^ + ^2 + b^cx 

+ b 5 aSf 

C m =c 0 + c x S e + c 2 <x 2 + c 3 S f + c 4 aS e 
+ c 5 aSj + c 6 a 
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Figure 1 Free-flying Airplane for Sub-scale 
Experimental Research (FASER) in the NASA Langley 
12 ft Low-Speed Wind Tunnel 
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Figure 2 Static Engine Thrust, Wind Tunnel Air Off 
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Figure 3 Angle of Attack Sweep for Subspace 
Boundary Definition, /?=0, pwr= 0 
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Sideslip Angle (deg) 


Figure 4 Sideslip Angle Sweep for Subspace 
Boundary Definition, 0^80, /nvr=0 



Power Level (percent) 


Figure 5 Power Sweep for Subspace Boundary 
Definition, cx= 40, J3= 0 



a (deg) 


Figure 6 Longitudinal Subspaces, Power Off 



a (deg) 


Figure 7 Longitudinal Subspaces, Power On 



a (deg) 


Figure 8 Lateral/Directional Subspaces, Power Off 
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Figure 9 Experiment Design Projection onto Two 
Dimensions 
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Figure 10 Predicted Squared Error (PSE) Components 
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Figure 13 Modeling and Prediction Results for 
Pitching Moment Coefficient, Longitudinal 
Subspace 16 
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Figure 14 3D Plotting Tool using Pitching Moment CofTicient for Longitudinal Subspace 16 
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